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Abstract 

The effects of 17-p-estradiol in osteoblasts are primarily mediated by the nuclear transcription factors, estrogen receptor 
(ER)a and ERp. ERs function through three general modes of action: DNA-binding dependent through estrogen response 
elements (EREs; designated nuclear ERE signaling); nuclear signaling via protein-protein interactions to other transcription 
factors (nuclear non-ERE signaling); and extra-nuclear signaling (membrane-bound functions of ERs). Identification of the 
specific transcriptional signatures regulated by each of these modes of action should contribute to an enhanced 
understanding of estrogen signaling in osteoblasts. To achieve this goal, we utilized specific mutations of ERa that eliminate 
the ability of the receptor to signal through a specific mode of action. The non-classical ERa knock-in (NERKI) mutation is 
incapable of signaling through direct DNA binding to EREs and the nuclear only ERa (NOER) mutation eliminates all 
membrane-localized signaling. Comparison of the gene expression patterns elicited by these mutations with the wild-type 
ERa (WT) pattern provides mode-specific data concerning transcriptional regulation by ERa. We expressed these constructs 
in the ER-negative osteoblastic cell line hFOB {—/+ estrogen) and performed global RNA-sequencing. Using a series of pair- 
wise comparisons, we generated three lists of genes that were regulated either by the nuclear ERE-dependent, nuclear ERE- 
independent, or extra-nuclear actions of ERa. Pathway and gene ontology analyses revealed that genes regulated through 
the nuclear ERE and nuclear non-ERE pathways were largely involved in transcriptional regulation, whereas genes regulated 
through extra-nuclear mechanisms are involved in cytoplasmic signaling transduction pathways. We also intersected our 
data with genes linked to bone density and fractures from a recent genome-wide association study and found 25 of 72 
genes (35%) regulated by estrogen. These data provide a comprehensive list of genes and pathways targeted by these 
specific modes of ERa action and suggest that "mode-specific" ligands could be developed to modulate specific ERa 
functionality in bone. 
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introduction 

Estrogen (E) signaling is mediated by two receptors, estrogen 
receptor (ER)a and ER|3. Although important for many tissue 
systems, ER signaling is of particular interest with regards to bone 
biology, as declining E levels during the menopause lead to 
increased bone turnover and bone loss, as is observed in 
postmenopausal osteoporosis [1]. Although extensive research 
has been conducted to understand the role of ERs in bone through 
use of cellular and animal models, the understanding of gene 
regulation and cellular pathways regulated by E in osteoblasts is 
stiU incomplete. 

Mechanistically, ERs modulate gene expression using a number 
of distinct modes of action. The classical pathway involves direct 
binding to estrogen response elements (EREs) in the control 
regions of genes, whereas an alternative nuclear pathway involves 
protein-protein interactions which are ERE-independent. To 
facilitate study of these pathways, an ERa mutation [non-classical 
ERa knock-in (NERKI)] has been made that eliminates ERE but 
retains non-ERE signaling (both nuclear and extra-nuclear) [2], 
allowing investigators to study the effects of non-ERE mediated 



ERa signaling alone. Mice harboring this mutation display an 
osteoporotic phenotype, demonstrating that nuclear ERE-mediat- 
ed signaling is important in regulating bone metabolism [3,4,5]. 

It has become increasingly appreciated that the nuclear actions 
of ERs cannot completely explain all aspects of ER signaling. 
Compelling evidence demonstrates that estrogen can act rapidly 
(within 5-15 minutes of E treatment) through non-genomic 
mechanisms, to regulate signal transduction through pathways 
such as ERK and PI3K [6,7,8,9]. This mechanism involves 
tethering of ERa to the plasma membrane (PM) via palmitoylation 
of a cysteine residue in its E domain, which facilitates association 
with caveolin- 1 [7] . A mutation of this cysteine has been made 
(termed "nuclear-only ER " or NOER) that completely eliminates 
membrane localization of ERa [10], and thereby prevents PM 
ERa signaling. Therefore, in comparison with wild-type ERa 
(WT), the NOER receptor can be used as a tool to identify those 
genes regulated through either the nuclear or membrane- 
associated ERa pathways. 

In an effort to identify target genes and pathways regulated by 
each of these distinct signaling mechanisms, we expressed WT, 
NERKI, or NOER receptors in an ER-negative osteoblast cell 
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model and measured global gene expression patterns following 
estrogen treatment using RNAseq. Comparison of estrogen- 
dependent gene expression patterns allowed us to compartmen- 
talize these patterns according to the known cellular mechanisms 
used by ERa, such as dependence on an ERE or whether these 
genes are regulated through the extra-nuck'ar (e.g. m(;mbrane- 
associated) function of ERa. These data identify characteristic E- 
regulated gene signatures from three distinct ERa regulatory 
mechanisms, which could potentially be used for targeting specific 
cellular ERa pathways. 

Materials and Methods 

Cell culture reagents, adenoviruses and infection of hFOB 
1.19 cells 

The hFOB 1.19 human fetal osteoblastic cell line (hFOB), 
produced by Harris and colleagues [11], was passaged in phenol 
red-free aMEM growth medium (Invitrogen, Carlsbad, CA) 
supplemented with IX antibiotic/antimycotic (Invitrogen), 10% 
(v/v) feted bovine serum (Hyclone, Logan, UT) and 300 ug/mL 
G418 selection antibiotic (Invitrogen). The Ad-ERa (expressing 
WT mouse ERa) and Ad-NERKI (expressing mouse ERa 
containing the double mutation E207A/G208A, which is incapa- 
ble of binding DNA [2]) adenoviral vectors were produced as 
previously described [12]. To produce the nuclear-only ERa 
(NOER), which exhibits solely nuclear signaling [10], the C451A 
mutation was introduced into WT ERa using the QuikChange II 
Site-Directed Mutagenesis Kit (Agilent Technologies, Santa Clara, 
CA) and subsequentiy used to produce an adenoviral vector 
(Vector Biolabs, Philadelphia, PA), resulting in "Ad-NOER". All 
ERa molecules were FLAG-epitope tagged at the N-terminus. 
hFOB cells were plated in 10-cm culture dishes (n = 6) at a final 
density of 2x10* cells/cm^. Cells were infected at the time of 
plating at a multiplicity of infection (MOI) of 15 (Ad-ERa), 30 (Ad- 
NERKI) and 22.5 (Ad-NOER), which is herein shown to resuk in 
approximately equivalent protein levels for ERa, NERKI and 
NOER (Fig. 1). The cells were infected in the presence of 8 |a,g/ 
mL hexadimethrine bromide (Sigma-Aldrich, St. Louis, MO) to 
enhance adenoviral infection. Following another 24 hour incuba- 
tion to allow for complete infection, the cells were treated with 
either ethanol vehicle (0.1% v/v) or 10 nM 17-fi-estradiol (0.1% 
v/ \'; Sigma-Aldrich) in the presence of media containing triple- 
stripped charcoal-treated FBS for an additional 24 hrs. 

Western Blotting 

Whole cell extracts from adenoviral-infected hFOB cells was 
prepared and subjected to western blotting analysis as previously 
described [12]. A total of 10 |ig protein, as determined by the 
Pierce BCA Protein Assay Kit (Thermo Scientific, Rockford, IE) 
was run from each adenoviral condition. The blots were incubated 
with primary antibodies directed against mouse-anti-Flag (1:1000) 
and mouse-anti-fi-actin (1:10,000). The mouse anti-IgG secondary- 
antibody hnked to horseradish peroxidase was used at 1:5000 
dilution. All antibodies were purchased from Sigma-Aldrich. Blots 
were visuahzed using enhanced chemUuminscence (GE Life- 
sciences, Piscataway, NJ) and exposed to X-ray film. 

RNA/cDNA isolation 

Total RNA was prepared using RNeasy minicolumns (Qiagen, 
Valencia, CA) and treated with RNase-free DNase (Qiagen) to 
remove potential contaminating DNA, as previously described 
[12]. The resulting RNA was used for either RNAseq analysis or 
in semi-quantitative real-time PCR analysis (QPCR). For the latter 
analysis, one microgram of total RNA was used in a reverse 



transcriptase (RT) reaction using the High Capacity cDNA 
Reverse Transcription Kit (Apphed Biosystems by Life Technol- 
ogies, Foster City, CA). The cDNA was diluted 1:5 with water 
prior to QPCR analysis. 

RNA sequencing (RNAseq) 

RNA libraries for RNAseq analysis were prepared from 1 00 ng 
of isolated total RNA from each sample using the manufacturer's 
instructions. Unique indexes were incorporated at the adaptor 
ligation step for loading multiple samples per flow cell. Three 
distinct indexed libraries were loaded per flow ceU and sequenced 
on an lUumina HiSeq 2000 using TruSeq SBS sequencing 
software (version 3) and SCS data collection software (version 
1.4.8). Base calling was performed using lUumina RTA (version 
1.12.4.2). An average of 130 million reads per sample was 
achieved resulting in ~96% mapped reads. A more detailed 
description of the RNAseq methodology used is as previously 
described [13]. The RNAseq data discussed in this publication 
have been deposited in NCBFs Gene Expression Omnibus [14] 
and are accessible through GEO Series accession number 
GSE55769 (http:/ /www.ncbi.nlm.nih.gov/geo/query/acc. 

cgi?acc = GSE55769). 

Semi-quantitative Real-time PCR Analysis (QPCR) 

The PCR reactions were run in the ABI Prism 7900HT Real- 
Time System (Applied Biosystems, Carlsbad, CA) using the 
Quantitect SYBR Green Master Mix (Qiagen), as previously 
described [15]. The method for data normalization using multiple 
reference genes and threshold calculations are as previously 
described [15]. Primer sequences for individual genes were 
designed using the Primer Express program (Applied Biosystems) 
and are available on request. 

Statistical analyses 

For the QPCR analyses, calculations and statistical comparisons 
were performed using Microsoft Office Excel 2003 (Microsoft 
Corp., Redmond, WA) and the data are presented as the mean ± 
SEM. All values of p<0.05 were considered statistically significant 
using two-sided t tests. The RNAseq data was analyzed essentially 
as previously described [13]. Briefly, paired-end reads from the 
raw RNAseq data were aligned using TopHat (version 2.0.6) [16] 
against the hl9 genome build using the bowtie 1 option [17], and 
quality control assessments were made using RSeQC software 
[18]. Gene counts were generated using HTSeq software and gene 
annotation files were obtained from lUumina. Gene expression 
data exhibiting a p-value^0.05, a false discovery rate of q£0.05 
and median gene counts in at least one group of 1 0, were used for 
further investigation. In order to identify which genes are 
regulated via non-ERE or extra-nuclear actions of ERa, the WT 
ERa data was compared the NERKI and/or NOER RNAseq 
datasets. Overlapping genes were identified using the "Venny" 
web resource (http:/ /bioinfogp.cnb. csic.es/tools/venny/). For the 
genes regulated via each of the three modes of action, pathway 
analysis was performed using the Ingenuity Pathway Analysis 
software (Ingenuity Systems, Redwood City, CA). Gene ontology 
(GO) terms from each dataset were generated using DAVID 
Bioinformatics Resources Version 6.7 [19]. 

Results 

RNAseq analysis of hFOB human osteoblastic cells 
expressing ERa modal variants 

The goal of this study was to identify and characterize estrogen- 
dependent gene expression patterns elicited by ERa through the 
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nuclear ERE-dependent, nuclear ERE-independent, or through 
extxa-nuclear signaling (e.g. membrane) in the ER-negative, 
human osteoblastic cell line hFOB [1 1] using global RNA 
sequencing (RNAseq). As tools to facilitate this approach, we 
utilized an ERa mutation that eliminates DNA binding through 
EREs (NERKI) [2,20] and an ERa mutation which can only 
signal through the nucleus (NOER) [10]. By comparing the gene 
expression patterns of these mutant ERa receptor with wild-type 
ERa (VVT), we can identify those patterns elicited by each mode of 
ERa action. Importantly, our analysis required that estrogen- 
regulated genes in each of the three modes of action were also 
regulated in WT, thus ensuring that these are physiologically 
regulated genes by ERa. Therefore, hFOB cells were infected with 
adenoviruses expressing FLAG-tagged WT, NERKI or NOER. 
Western blot analysis confirmed that all ERs were expressed at 
~equal levels (Figs. lA/B). 

To identify differentially expressed genes with each of the ERa 
modal variants, hFOB cells were infected with adenoviral vectors 
expressing WT, NERKI or NOER receptors, treated with 
estrogen (10 nM E2) for 24 h and RNAseq was performed (see 
Materials and Methods for details). In a previous study, we had 
determined that genes exhibiting a median gene count of at least 
10 in one group represent an "expressed" gene [13]. Therefore 
genes with a median gene count of less than 10 in both comparison 
groups were called "non-expressed" and excluded from further 
analysis. In this study, the remaining genes exhibiting a p-valueS 
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Figure 1 . Estrogen receptor expression in liFOB cells. A) hFOB 
cells were infected with Ad-ERa, Ad-NERKI and Ad-NOER and cultured 
for 24 h. Protein extracts were prepared and a western blot was 
performed using the anti-FLAG and anti-P-actin antibodies. B) 
Densitometry was performed and the data are expressed as ER 
expression relative to the p-actin control. 
doi:1 0.1 371/journal.pone.0095987.g001 



0.05 and a false discovery rate of qS0.05 between control and 
estrogen treatments were used for further analyses. 

Using these criteria, an estrogen-regulated gene Kst for WT was 
generated. In the WT treatment groups where all modes of 
estrogen action are preser\'ed, a total of 4353 estrogen-regulated 
genes were identified (Table SI). Thirty randomly chosen genes 
from this dataset were analyzed by Q_PCR to validate the RNAseq 
data and a high correlation (r = 0.98, P<0.001) was observed 
(Figure SI), demonstrating a high level of accuracy of the RNAseq 
technology. 

The flow-chart depicted in Figure 2 describes our data analysis 
strategy in detaU. SimUarly to WT, the NERKI and NOER 
datasets were passed through our data filters (median gene counts 
10, p<0.05 and q<0.05) and 2375 estrogen-regulated genes were 
identified from the NERKI dataset (termed "ERE-independent") 
and 3147 estrogen-regulated genes were identified from the 
NOER dataset (termed "Nuclear-only) (Fig. 2A). Complete lists for 
all three datasets can be found in Table SI (tabs 1-3, respectively). 
To identify the "ERE-dependent" genes, the 2375 genes from the 
NERKI dataset were subtracted from the 4353 genes from the 
WT dataset (Fig. 2B; top row). This subtraction results in 1978 
genes that are regulated by WT but not regulated by NERKI, 
therefore by inference these genes are regulated via ERE- 
dependent mechanisms. Similarly, to identify those genes whose 
estrogen-dependent regulation occurs from signals generated 
outside the nucleus (termed "Extra-nuclear), the 3147 genes in 
the Nuclear-only dataset were subtracted from the 4353 genes 
from the WT dataset (Fig. 2B, bottom row). This subtraction 
results in 1206 genes that are regulated by WT, but not the 
NOER, suggesting these genes are regulated in part by 
mechanisms originating outside of the nucleus. However, our 
ERE-dependent and ERE-independent gene lists include regula- 
tion occurring in both the nucleus and outside of the nucleus. To 
identify nuclear genes that are ERE-dependent and ERE- 
independent, each of these lists were overlaid with the 1 206 genes 
in the Extra-nuclear list. For the ERE-dependent genes, this 
resulted in 1167 genes termed "Nuclear ERE-dependent" and 
with 811 genes from the Extra-nuclear list (Fig. 2C; left side). 
Similarly, for the ERE-independent genes, this resulted in 1980 
genes termed "Nuclear ERE-independent" with 395 genes from 
the Extra-nuclear list (Fig. 2C; right side). Note that the sum of the 
Extra-nuclear genes from both comparisons is the complete Extra- 
nuclear gene list (81 1+395 = 1206). In this way, the lists comprising 
the three modes of ERa action (Nuclear ERE-dependent, Nuclear 
ERE-independent, Extra-nuclear) were identified. Complete lists 
for these datasets can be found in Table S2 (tabs 1-3, respectively). 
The percentage breakdown and number of estrogen regulated 
genes from each dataset representing these distinct modes of ERa 
action is shown in Figure 3. The top 20 up- and down-regulated 
genes for each dataset are shown in Tables 1-3. 

Ingenuity pathway analysis (IPA) and gene ontology (GO) 
analysis of the datasets from the three modes of ERa 
action 

Following global identification of estrogen-regulated genes 
comprising the three modes of ERa action (described above), we 
performed IPA on each dataset to identify those cellular pathways 
mostiy like regulated in each of these modes (see Table S3 for 
complete lists). We focused on those pathways with known 
regulatory roles in osteoblast biology. Pathways with known roles 
in bone biology' for the Nuclear ERE-dependent (Table 4), 
Nuclear ERE-independent (Table 5) and Extra-nuclear (mem- 
brane) (Table 6) are indicated. 
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Figure 2. Flow-chart detailing our data analysis strategy (described in detail in the text). 

doi:1 0.1 371 /journal.pone.0095987.g002 




Figure 3. Confirmation of the RNAseq datasets and percentage 
breakdown of all estrogen-regulated genes among the three 
modes of ERa action. Breakdown of all the estrogen-regulated genes 
in the RNAseq analysis among the modes of ERa action, showing 
percentages and the number of genes from each mode in parentheses. 
doi:1 0.1 371 /journal.pone.0095987.g003 



We also performed gene ontology (GO) analysis to categorize 
the target genes regulated in each mode of ERa action into specific 
molecular functions (GOTERM_MF_FAT) or cellular compo- 
nents (GOTERM_CC_FAT) (see Table S4). We found diat the 
target genes regulated in the nuclear ERE and non-ERE mode of 
ERot action were largely enriched for genes involved in 
transcriptional regulation. Interestingly, the target genes regulated 
by the extra-nuclear mode of ERot action, presumably occurring 
through membrane-associated processes, were largely enriched for 
genes involved in metabolic pathways (i.e. NADH/NADPH 
oxidoreductase activity and FAD binding) and structural compo- 
nents of the ribosome. 

Comparison of the datasets representing the three 
modes of ERa action with genes with known 
involvement in bone biology 

To place our results in a more physiologically relevant context 
in terms of bone biology, we intersected our gene lists 
corresponding to the three modes of ERot action with a recent 
meta-analysis that identified 72 genomic loci with genome- wide 
significance to bone mineral density and/ or fracture rates [2 1] . As 
seen in Figwe 4, a total of 25 (—35%) of the GWAS-identified 
genes are found in our datasets. 
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Table 1. Most Highly Regulated Genes in the Nuclear ERE-dependent Dataset. 



Up-regulated Genes 






Down-regulated Genes 






Gene Symbol 


WT 


NOER 


Gene Symbol 


WT 


NOER 


PGLYRP2 


103.00 


5.56 


GALNT13 


0.22 


0.32 


CI lorfse 


25.39 


13.59 


SLITRKl 


0.26 


0.30 


MIR 1286 


17.33 


122.68 


KCNJ2 


0.27 


0.41 


TH 


15.79 


11.60 


C6orfl41 


0.33 


0.41 


RNF223 


15.72 


7.38 


DUSP6 


0.36 


0.33 


TMIE 


12.46 


36.10 


TNF 


0.36 


0.47 


NKX2-3 


10.27 


2.72 


UNC00158 


0.36 


0.31 


AQP5 


8.70 


8.02 


CLDN14 


0.38 


0.50 


IGL0N5 


7.35 


5.04 


C18orfl 


0.39 


0.29 


OXT 


6.90 


9.97 


KLF14 


0.42 


0.47 


WISP2 


6.50 


8.57 


MIR125B1 


0.43 


0.23 


MYTIL 


6.37 


7.11 


VIT 


0.43 


0.63 


HS3ST6 


6.27 


4.00 


LOCI 00506895 


0.44 


0.51 


GFAP 


6.26 


3.52 


IL6 


0.45 


0.34 


FOXES 


6.20 


4.25 


GMNC 


0.48 


0.57 


MlOX 


6.13 


3.79 


E5M1 


0.50 


0.57 


CDX2 


6.08 


12.87 


SLITRK6 


0.51 


0.47 


FGF8 


6.01 


7.18 


ARL14 


0.51 


0.55 


LOC389332 


5.08 


7.58 


RUNXni 


0.51 


0.70 


LOC728228 


4.93 


3.97 


MYCTl 


0.52 


0.45 



doi:l 0.1 371 /journal.pone.0095987.t001 



Discussion 

Estrogen receptor (ER)-a can utilize a number of difTerent 
mechanisms to modulate gene expression, including nuclear ERE- 
dependent, nuclear ERE-independent, and extra-nuclear signaling 
(e.g. membrane). The goal of this study was to identify the global 
gene expression patterns regulated by each of these mechanisms in 
the ER-negative, human osteoblastic cell line hFOB [1 1] using 
global RNA sequencing (RNAsecj), in an effort to further 
understand how each mode of ERa action contributes to the 
overall response of osteoblasts to estrogen. Using specific ERa 
mutants which eliminate one mode of ERot action, and by 
comparing these datasets to the WT dataset, we were successful in 
identifying genes regulated by either the nuclear ERE-dependent, 
nuclear ERE-independent, and extra-nuclear actions of ERot. 

Examination of the relative proportions of the modes of ERa 
action revealed that a majority of estrogen-regulated genes (72%) 
are regulated at least partially via nuclear mechanisms. This data is 
supported by recent report that 82% of estrogen-regulated genes 
in the uterus are regulated via nuclear ERa pathways [22]. Of 
these, nearly two-thirds (63%) are from ERE-independent 
mechanistic pathways and one-third (37%) from ERE-dependent 
mechanistic pathways. This suggests that a common mechanism 
for ERa to influence estrogen-dependent transcription in the 
nucleus is through cooperation with other DNA-bound transcrip- 
tion factor complexes. Our pathway analysis of each mode of ERa 
action supports this contention, since 135 cellular pathways were 
significantiy regulated by the nuclear ERE-independent dataset, 
whereas only 48 and 33 pathways were sigiiificandy regulated in 
the nuclear ERE-dependent and extra-nuclear datasets, respec- 
tively. 



Recent research has demonstrated that ERa is found at the 
plasma membrane of nearly all ERa-expressing cells [10], and 
therefore contributes to overall estrogen signaling in most ER- 
positive cells, including osteoblasts. Indeed, we find that 28% of all 
estrogen-responsive genes originate from extra-nuclear (e.g. 
membrane) compartment, since they are regulated by WT ERa 
but not the nuclear-only ERa (e.g. NOER). In addition, 
accumulating evidence now suggests that significant cross-tall^ 
occurs between the membrane and nuclear ERa signaling 
[10,23,24,25]. A well-known example of this cross-talk is in the 
estrogen-dependent regulation of the Ccndl gene in breast cancer 
cells, where membrane-initiated activation of PI3K/AKT and 
ERK signaling, as well as ERE-dependent mechanisms are 
necessary for maximal estrogen-stimulated transcription [7,26,27]. 

An important and interesting aspect of this research is that 
estrogen-regulated pathways that are either in common or unique 
to the three modes of ERa action can be identified. Comparison of 
the pathways from all modes of ERa action revealed a few shared 
pathways, such as: integrin and integrin linked kinase (ILK) 
signaling, tight junction signaling, IL8 signaling, MAPK mTOR 
signaling. These may represent a more generalized function of 
ERa action, in which multiple modes are involved. However, most 
of the significant pathways are unique to each mode, suggesting a 
unique biology is involved dependent on how ERa is functioning. 
Unique pathways regulated in the nuclear ERE-dependent mode 
include EGF signaling, sonic hedgehog signaling and many 
cholesterol and nucleotide biosynthetic pathways. This is in 
contrast to the more common signaling pathways regulated m 
the nuclear ERE-independent mode of ERa action such as PPAR, 
BMP, Wnt, GR, IL6, TGFP, among many others. Unique 
signaling pathways in the extra-nuclear mode of ERa action 
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Table 3. Most Highly Regulated Genes in the Extra-nuclear Dataset. 



Up-regulated Genes 




Down-regulated Genes 




Gene Symbol 


WT 


Gene Symbol 


WT 


WNT3A 


84.99 


ADCY2 


0.07 


UNCX 


45.10 


MIR483 


0.12 


LEFTYl 


15.34 


SLC32A1 


0.15 


P0M121L9P 


10.56 


LINC00189 


0.17 


0T0P2 


10.33 


KRT84 


0.18 


KLK4 


7.42 


ZIC5 


0.23 


AP0A4 


6.94 


ACSM4 


0.26 


PSAPLl 


6.03 


MIRM69 


0.29 


MYCN 


5.65 


HIST1H3A 


0.30 


TLX! 


5.59 


TCEAL2 


0.30 


CPN2 


5.23 


LOC100129046 


0.31 


PPP1R27 


5.21 


KRT5 


0.35 


LCN2 


5.15 


ME0X2 


0.35 


MYH6 


4.95 


CPR88 


0.37 


INSMl 


4.91 


ISU 


0.37 


LOC727710 


4.81 


BUD 


0.39 


Clorfei 


4.73 


KCNJ2-AS1 


0.39 


MB 


4.70 


LOC401 164 


0.39 


CHRND 


4.40 


RNU6ATAC 


0.40 


LOC100507218 


4.33 


SPINLWl 


0.41 
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include FAK, AMPK, calcium and cNOS signaling. Understand- 
ing how each mode regulates a unique aspect of overall ERa 
signaling in osteoblasts may lead to the generation of "mode- 
specific" ligands used to target specific pathways in estrogen- 
responsive cells. 

At a physiological level, estrogen signaling is well recognized as 
an important determinant of bone mineral density (BMD), as 
declining estrogen levels during the menopause lead to increased 
bone turnover and bone loss, as is observed in postmenopausal 
osteoporosis [1]. Therefore it is not surprising that significant 
overlap exists upon intersection of our estrogen-dependent gene 
lists with the 72 loci identified in GWAS studies as being related to 



BMD and/or fracture rates [21]. These include genes with well- 
known roles in osteoblast biology, such as Rmx2, Sppl [osteopontiri), 
Esrl (ERa.) and Tnfr.sfllh (Opg), as well as other genes with less 
recognized roles in bone. This demonstrates that all modes of ERa 
action contribute to gene regulatory patterns important in bone 
biology. Future studies are needed to examine the role of these 
genes in bone diseases such as osteoporosis. 

In reality, these modes of ERa action do not exist iii isolation. It 
is well known that the rapid effects of membrane estrogen 
signaling have modulatory effects on nuclear estrogen signaling 
through phosphorylation of ERa and/ or coregulatory transcrip- 
tion factors, as well as induction of numerous signaling cascades 



Table 4. Selected Pathways of Interest with Known Function in Bone in the "Nuclear ERE-dependent" Dataset. 





Top Canonical 
Pathways 


/'Value 


Ratio 


Genes 


GNRH Signaling 


0.0081 


0.11 


MAP3K9,SRC,MAP3K6,MAPK8,DNM3,DNMl PAK3,CREB1,!TPR3,PRKAR1B,PRKCE,PRKCH,MAP2K3,ELK1,MAP2K1 


EGF Signaling 


0.0166 


0.14 


SRC,!TPR3,MAPK8,PIK3CD,CSNK2B,ELK1, MAP2K1,RASA 1 


IL-8 Signaling 


0.0240 


0.09 


RAC2,SRC,PLD2,MAPK8,HBEGF,UMK2,BAX, CSTB,EIF4EBP1,CCND3,PRKCE,PIK3CD, PRKCH,GNA13,KDR,RH0F,MAP2K1,IRAK2 


Sonic Hedgehog 
Signaling 


0.0263 


0.16 


ADRBKl,PTCHhPRKARlB,HH!P,DYRKlA 


VDR/RXR 
Activation 


0.0380 


0.12 


SPP1,CCNC,RUNX2,IURL1,MXDIPRKCE PRKCH,SEMA3B,CYP27B1 


PDGF Signaling 


0.0380 


0.11 


SRC,MAPK8,PDGFRA,PIK3CD,CSNK2B,ELK1, MAP2KIRASA IJNPPSD 


p38 MARK 
Signaling 


0.0427 


0.10 


TGFBR2,DAXKATF1,IURL 1,MEF2D,CREB 1, MEF2A,MAP2K3,ELK1,TNF,FAS,IRAK2 
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Table 5. Selected Pathways of Interest with Known Function In Bone In the "Nuclear ERE-independent" Dataset. 



Top Canonical Pathways 


/'Value 


Ratio 


Genes 


Role of Osteoblasts, 
Osteoclasts and Chondrocytes 
in Rheumatoid Arthritis 


0.0000 


0.22 


BMP4, NFATC3, PIK3R1, TCF7, BMP8A, TGFBl, PIK3CG, WNT78, WNT4, SMADl, ADAMTS5, BMP88, 
TNFRSniB, WNT9A, PPP3CC, MAPKU, PIK3R3, IU8, SFRP1, RELA, NFKBIE, FZDl, NFKBl, NFATCl, SMURFl, 
JUN, NFKBIA, BMPRIA, NGFR, TNFRSFIB, TNFRSFIA, DVLl, DLXS, IKBKE, NFATC4, TNFRSFl lA, CALMl, FZD8, 
FOB, TRAF2, FZD4, FOXOl, CSFl, ILIB, DKKl, BMP6, WNTl 1, LRPl, TCF7L2 


PPAR Signaling 


0.0000 


0.26 


PPARA,RELAPDGFANFK8IE,NFKBhNR2Fl, NFKBIA,JUN,NGFR,MRAS,TNFRSF18,aTED2, 
TABTJNFRSFUB.SRAl, TNFRSFl A.IKBKE, PDGFB,AIP,F0S,TRAF2,IL18,IL1B,NC0R2, PTGS2,RXRA 


Wnt/p-catenin Signaling 


0.0000 


0.22 


CDKN2A,TGFBR3,PPP2R5B,GSK3A,FZDh RARG,CCND1,TCF7,JUN,NLK,TGFBIWNT7B, 
RARA,TGFB2,WNT4,SOXW,TABl,SOX4,SOX7, CSNK1G2,WNT9A,DVL1,S0X1 hACVRlB,FZD8, 
FZD4,TLE4.TLE3,SOX8,DKKlPINlSFRPh DVL2.ACVR2A,WNT11,LRP1,TCF7L2 


IL-6 Signaling 


0.0000 


0.23 


RELA,S0CS3,S0CS 1,NFKBIE,PIK3R I.NFKB 1, VEGFA,JUN,NFKB!A,NGFR,P!K3CG,MRAS, 
TNFRSFl B.MAPKAPK2,TABhTNFRSFllB. MCL1,IL8,TNFRSF1A,IKBKE,CEBPB,MAPK12, 
PIK3R3,F0S,1L18,TRAF2,IL 1B,HSPB7 


TGF-p Signaling 


0.0000 


0.25 


RUNX3,BMP4,SMAD3,SK!,SMAD7,PITX2, MAPK12,ACVR1B,!NHB8,!NHBA,SMURF1,F0S, 

ni AC /t II IKI DA Ann 1 A 'T/'~CD 1 A An AC T/~nO 1 TCCD C IIA A ^/"l/DT A TA O f 

ri/\b4,JUN,DMrn lA, 1 bro l,MHAt>, 1 \jro2, 1 rtJibMAu l,ALvH2A,IAtS I 


BMP signaling pathway 


0.0003 


0.23 


RELA,BMP4,FST.SMAD7,NFKBIMAPK12, PITX2,N0G,SMURFIJUN,BMP8A,BMPR1A, 
MHAS,PHl\AG2,DlVlro,SIVtAU l,DmH8D, lABl 


VDR/RXR Activation 


0.0005 


0.23 


IL12A,PDGFA,HES1,CEBPB,TH8D,HR,KLF4, GTF2B,GADD45A,FOX01,PRKCD,HOXA10, 
CDKN1A,IGFBP3,CEBPA,TGFB2,NC0R2,RXRA 


NF-kB Signaling 


0.0014 


0.17 


RELA,BMP4,NFKB!EJGFBR3,P!K3R1, NFKBl, NGF,FGFR3,NFK8IA,CARD10,8MPR1A,NGF, 
PIK3CG,NTRK1,MRAS,TNFRSF1B,TAB1, TNFRSFl 1B,TNFRSF1A,RELB,MALT1, 
TNFRSFl 1A,TLR9,PIK3R3,TLR4,IL18,TRAF2, TGFA.ILIB 


Notch Signaling 


0.0028 


0.24 


DLL 1,JAG2,MFNG,LFNG,HES IJAG 1,DLL4, NOTCH 1,DTX2,HEY1 


IGF-1 Signaling 


0.0032 


0.19 


S0CS1,S0CS3,IGFSP4,CTGF,JAK1,PIK3R1, SOCS2,SOCS6,PIK3R3,FOS,JUN,FOXOl 
PIK3CGJRS 1,IGFBP3,MRAS,PRKAG2,IRS2, CYR6 1 


Estrogen Receptor Signaling 


0.0457 


0.14 


MED12L,PELP1,SRA1,TAF6,PCK1,TAF10, ERCC2MED12,G6PC3,TAF6L,GTF2B,POLR2A, 
MED17,GTF2H4,GTF2El,MRAS,MED16,NCOR2, ESR 1 
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(reviewed in [6,28]). In this way, elimination or modulation of any 
one mode of ERa action affects the whole system, eventually 
affecting gene expression and phenotypic expression of the cell, 
and ultimately the organ and organism [9] . These actions may be 
very different in various estrogen-responsive tissues, depending on 
the amount of ERa present and the complement of transcription 
factors expressed in that particular cell. 

Since we treated the ER-infected hFOB cells for 24 h, our data 
most certainly contains not only genes that directly respond to 
ligand-bound ERa (a primary response), but also those genes 
whose regulation is dependent on primary response genes (e.g. 



secondary response genes). Although this is an important aspect to 
consider when examining the gene lists produced in this study, our 
primary premise is that these 3 modes of ERa action which are 
mechanistically distinct, ultimately put into motion a mode- 
specific program that leads to the regulation of distinct subsets of 
genes. Future studies with shorter treatment periods will be 
necessary to limit these lists to primary response genes. 

We recognize that our study has several limitations. First, we are 
using an in vitro cell culture system and since the hFOB cells are 
ER-negative, we expressed these ER mutations using an adeno- 
viral vector. Using an ER-negative cell model was essential to 



Table 6. Selected Pathways of Interest with Known Function in Bone in the "Extra-nuclear" Dataset. 





Top Canonical 
Pathways 


Value 


Ratio 


Genes 


Integrin Signaling 


0.0014 


0.11 


RAP1B,PIK3CAJSPAN5,MPRIP,DIRAS3,ACTA2,A8L1,P!K3R5,TSPAN2,RH0J,PIK3R4,PTEN, 
MYL9,ITGAE,ITGA3,PAKhCAPNSl,GRB7 RHOU,ACTC1,ITGB5,CAPN10 


mTOR Signaling 


0.0014 


0.11 


NAPEPLD,PIK3CA,PLD3,PRKAB1,DIRAS3, PIK3R5,PPP2R3B,RPSS,VEGFS,RHOJ,P!K3R4, 
RICTOR,EIF4Gl,RPS7,PRKAA2,EIF3A,RHOU, PRRS.PPP2R2C.INSR.RPS14 


FAK Signaling 


0.0087 


0.11 


PIK3CA,ITGA3,PAK1,CAPNS1,ACTA2,P!K3R5, PIK3R4,TNS1,ACTC1,CAPN10,PTEN 


AMPK Signaling 


0.0148 


0.10 


PIK3CA,PRKAS1,PIK3R5,PPP2R3S,LIPE,P!K3R,ADRB1,FASN,PRKAA2,PPP2R2C,CPT1C, MLYCD,INSR,CHRNA3 


Calcium Signaling 


0.0155 


0.09 


RAP18,MYH6,ACTA2,HDAC10,TPM3,MYH7, CHRND,MYH7B,GRIA4,MYL9,CAMK2A,MYH2, 
CASQ1,ASPH,ACTC1,CHRNA3,CAMK2G 


p53 Signaling 


0.0363 


0.11 


PIK3CA,THBS l,E2Fl,PIK3R5,C12orf5.MDM2, PIDD,PIK3R4,PTENJP53I3 


eNOS Signaling 


0.0468 


0.09 


8DKR82,NOSIP,ADCY2,PIK3CA,CCNA 1, HSPA 14,PRKAB1,PRKAA2,PIK3R5,VEGFB, PIK3R4,CHRNA3 


TR/RXR Activation 


0.0490 


0.10 


SLC16A3,PIK3CA,ADRB1,C0L6A3,FASN, PIK3R5,NCOA4,MDM2,PIK3R4 
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Figure 4. Comparison of tKie RNAseq datasets to gene witli l<nown involvement in bone biology from genome-wide association 
studies (GWAS). The nuclear ERE-dependent, nuclear ERE-independent, and extra-nuclear datasets were intersected with the 72 genes identified by 
GWAS studies. The overlapped genes for each mode of ERoc action are listed and color-coded for clarity. The directionality of the estrogen-dependent 
regulation is denoted with an arrow following the gene symbol. 
doi:1 0.1 371 /journal.pone.0095987.g004 



eliminate any endogenous ER expression that could confound the 
results. In addition, during development of this cell culture system 
we carefully titrated the various adenoviral vectors and used the 
minimal amount to achieve equal expression of each receptor in > 
95% of the cells. Further, while we have not formally compared 
the levels of ERa expression in our transfected cells with a panel of 
WT osteoblasts, it is likely that these levels are considerably higher 
in our transfected cells, despite the fact that we tried to use fairly 
low viral titers for the transfections. This issue is compounded by 
the fact that levels of ERot expression in vitro vary considerably with 
the stage of osteoblast differentiation. Thus, while the strength of 
our in vitro system is that we were able to study the various estrogen 
signaling pathways in isolation, we recognize that further 
validation of these findings in in vivo systems with physiological 
levels of ERa are needed. A second limitation is that by using 
hFOB cells, we limited our study to osteoblasts. However, these 
datasets may be useful in identifying similar ERa action in other 
estrogen-responsive tissues, such as breast and uterus. Further 
studies are needed to examine the genes/pathways regulated by 
these modes of ERot action in other estrogen-responsive tissues. A 
third limitation is that only ERot, and not ERP was studied. These 
modes of action also exist for ER(3, including the membrane- 
associated ERP mechanistic pathway that has been extensively 
studied in cardiomyocytes [29,30]. Future studies wiU be needed to 
clarify these modes of ERP action and interaction with ERot, 
which win be considerably more complex. 

In conclusion, we have identified three molecular modes of ERa 
action in human osteoblastic cells using RNAseq methodology. 
The identification of estrogen-dependent gene expression patterns 
and pathways within these modes of action is an important first 
step into understanding the complete picture of how ERa action 
contributes to overall ERot signaling in bone, and may provide the 
data necessary for generation of "mode-specific" ligands to 
preferentially influence a single mode of ERa action. 



Supporting Information 

Figure SI QPGR confirmation of the RNAseq dataset 
from wild-type ERa. Thirty randomly chosen genes from the 
wild-type ERa dataset were chosen and QPCR was performed on 
independent samples. The plot shows a high degree of 
concordance between the RNAseq dataset and the QPCR 
analysis, with an r = 0.98 and P<0.001. 
(TIF) 

Table SI Complete gene lists of the estrogen-regulated genes 
from the 3 ERa variants in this study. The estrogen-regulated 
genes from the wUdtype ERa, ERE-dependent, and Nuclear-only 
ERa variants (from Fig. 2A) are listed in each respective tab. The 
gene name, accession numbers, fold-change with estrogen (EC), p- 
value and false discovery rate (FDR) are indicated. 
(XLSX) 

Table S2 Complete gene fists of the estrogen-regulated genes 
comprising the 3 modes of ERa action. The "Nuclear ERE- 
dependent", "Nuclear ERE-independent" and "Extra-nuclear" 
estrogen regulated genes are indicated (from Fig. 2C) with the 
same values as described in Table SI. 
(XLSX) 

Table S3 Complete Ingenuity pathway analysis (IPA) of the 
genes from the 3 modes of ERa action. The IPA pathways for the 
genes for each mode of ERa action are indicated. 
(XLSX) 

Table S4 Gene ontology (GO) analysis of the genes from the 3 
modes of ERa action. Gene ontology (GO) analysis was performed 
to categorize the target genes regulated in each mode of ERa 
action into specific molecular functions (GOTERM_MF_FAT) or 
ceUular components (GOTERM_CC_FAT). 
(XLSX) 
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